#include "TStopwatch.h"

#include "AnaND280.h"

#include "load.h"

#include "bunchtime.h"

#include "NeutrinoUtils.h"


TChain * AnaND280::freconTracker     =0x0;
TChain * AnaND280::freconGlobal      =0x0;
//TChain * AnaND280::fBasicHeader      =0x0;
//TChain * AnaND280::fBasicDataQuality =0x0;
//TChain * AnaND280::fBeamSummaryData  =0x0;

TChain * AnaND280::fmcTruthRooTrackerVtx =0x0;
TChain * AnaND280::fmcTruthVertices      =0x0;
TChain * AnaND280::fmcTruthTrajects      =0x0;

//Int_t AnaND280::fEventID=-999;
Int_t AnaND280::fNTVtx=-999;
TClonesArray *AnaND280::fVertices =0x0;
Int_t AnaND280::fNTraj=-999;
TClonesArray *AnaND280::fTrajectories =0x0;
Int_t AnaND280::fNVtx=-999;
TClonesArray * AnaND280::fVtx =0x0;
Int_t AnaND280::fNPIDs=-999;
TClonesArray * AnaND280::fPIDs =0x0;
Int_t AnaND280::fNTracks=-999;
TClonesArray * AnaND280::fTracks =0x0;

Int_t AnaND280::ftTidPassMu=-999;
Int_t AnaND280::ftRidPassMu=-999;

TList * AnaND280::fList = 0x0;
//============================================================================

void AnaND280::SetNeutrinoRec()
{
  delete NeutrinoUtils::fNeutrinoRec;

  //beam direction - nd280, assuming along z axis
  NeutrinoUtils::fNeutrinoRec = new TVector3(0,0,1);
}

Bool_t AnaND280::FGDCut(const TLorentzVector * fropos)
{
  const Double_t FGD1xmin = -874.51;
  const Double_t FGD1ymin = -819.51;
  const Double_t FGD1zmin = 136.875;
  const Double_t FGD1xmax = 874.51;
  const Double_t FGD1ymax = 929.51;
  const Double_t FGD1zmax = 446.955;
  
  if(FGD1xmin < fropos->X() && 
     fropos->X() < FGD1xmax && 
     FGD1ymin < fropos->Y() && 
     fropos->Y() < FGD1ymax && 
     FGD1zmin < fropos->Z() && 
     fropos->Z() < FGD1zmax){
    return kTRUE;
  }
  else{
    return kFALSE;
  }
}

void AnaND280::IsTPCOK(const ND::TGlobalReconModule::TGlobalPID *trTrack, Bool_t &istpcok, Bool_t &ismuon)
{
  const ND::TGlobalReconModule::TTPCObject *tmptpc = 0x0;
  Int_t nfound = 0;

  for(Int_t itpc = 0; itpc < trTrack->NTPCs; itpc++){
    const ND::TGlobalReconModule::TTPCObject *gtpctrack =(ND::TGlobalReconModule::TTPCObject*)(*(trTrack->TPC))[itpc];
    //get values when in tpc2 for mu sel
    if(gtpctrack->Detector == 2){
      nfound++;
      tmptpc = gtpctrack;
    }
  }
  
  if(nfound>1){
    printf("AnaND280::IsTPCOK nfound>1 %d\n", nfound); exit(1);
  }

  istpcok = nfound;

  //--
  if(nfound){
    const Double_t charge = tmptpc->Charge;
    const Double_t pullm  = tmptpc->PullMuon;
    const Double_t pulle  = tmptpc->PullEle;
    const Int_t nodes  = tmptpc->NNodes;
    
    const bool passchrg = charge < 0;
    const bool passmuonpid = pullm < 2.0 && pullm > -2.0 && pulle < -2.0;
    const bool passtrqual = nodes > 18;
    
    ismuon = passchrg && passmuonpid && passtrqual;
  }
  else{
    ismuon = kFALSE;
  }
}

void AnaND280::SetNeutrinoMuonSim()
{
  //test recheck: algo may be problematic!!

  Int_t nfound = 0;

  const ND::TTruthVerticesModule::TTruthVertex *  nuout = 0x0;
  const ND::TTruthTrajectoriesModule::TTruthTrajectory * muonout = 0x0;
  for(Int_t ivt = 0; ivt < fNTVtx; ivt++){
    const ND::TTruthVerticesModule::TTruthVertex *vtxTrue = (ND::TTruthVerticesModule::TTruthVertex*)(*fVertices)[ivt];
    if(vtxTrue->ID != GettTidPassMu()){
      continue;
    }
    //get the PDG code for the track - have to match track true ID with a trajectory ID
    for(Int_t itr = 0; itr < fNTraj; itr++){
      const ND::TTruthTrajectoriesModule::TTruthTrajectory *trjTrue = (ND::TTruthTrajectoriesModule::TTruthTrajectory*)(*fTrajectories)[itr];
                        
      if(trjTrue->PDG == NeutrinoUtils::kPDGMuonMinus && trjTrue->ID == GettRidPassMu()){
        nuout = vtxTrue;
        muonout = trjTrue;
        nfound++;
      }
    }         
  }

  if(nfound>1){
    printf("AnaND280::SetNeutrinoMuonSim nfound>1!! %d\n", nfound); exit(1);
  }

  delete NeutrinoUtils::fNeutrinoSim;
  NeutrinoUtils::fNeutrinoSim = 0x0;

  delete NeutrinoUtils::fMuonSim;
  NeutrinoUtils::fMuonSim = 0x0;

  if(nfound==1){
    //confirmed with the one from  ND::NRooTrackerVtx *rooVertex 
    NeutrinoUtils::fNeutrinoSim = new TLorentzVector(nuout->NeutrinoMomentum * 1e-3);

    NeutrinoUtils::fMuonSim = new TLorentzVector(muonout->InitMomentum * 1e-3);

    //test
    //printf("test21 %f %f %f\n\n",NeutrinoUtils::fNeutrinoSim->E(), NeutrinoUtils::fNeutrinoSim->Vect().Phi()*TMath::RadToDeg(), NeutrinoUtils::fNeutrinoSim->Vect().Theta()*TMath::RadToDeg());

    if(nuout->NeutrinoPDG != NeutrinoUtils::kPDGNumu){
      printf("AnaND280::SetNeutrinoMuonSim warning nuout->PDG != NeutrinoUtils::kPDGNumu %d %d\n", nuout->NeutrinoPDG, NeutrinoUtils::kPDGNumu);
    }
  }

  NeutrinoUtils::SetMuonNSim(nfound);
}

void AnaND280::SetNeutrinoParentPDG(const Int_t ib)
{
  Int_t nfound = 0;
  NeutrinoUtils::NeutrinoParent parid = NeutrinoUtils::kOther;

  for(Int_t ivt = 0; ivt < fNVtx; ivt++){
    const ND::NRooTrackerVtx *rooVertex = (ND::NRooTrackerVtx*)fVtx->At(ivt);
    if(ib != BunchTime::GetVTBunchNo(rooVertex->TimeInSpill)){
      continue;
    }
  
    if(GettTidPassMu() == rooVertex->TruthVertexID){
      nfound++;

      const Int_t beampdg = rooVertex->NuParentPdg;
 
      if(beampdg == NeutrinoUtils::kPDGPionPlus || beampdg == NeutrinoUtils::kPDGPionMinus){
        parid = NeutrinoUtils::kPion;
      }
      else if(beampdg == NeutrinoUtils::kPDGKaonPlus || beampdg == NeutrinoUtils::kPDGKaonMinus || beampdg == NeutrinoUtils::kPDGKLong){
        parid = NeutrinoUtils::kKaon;
      }
      else if(beampdg == NeutrinoUtils::kPDGMuonMinus || beampdg == NeutrinoUtils::kPDGMuonPlus){
        parid = NeutrinoUtils::kMuon;
      }

      /*
      //test
http://www.hep.lancs.ac.uk/nd280Doc/oaAnalysisFormat/5E/ND__NRooTrackerVtx.html
        intStdHepPdg[100]pdg codes (& generator specific codes for pseudoparticles)
        intNEipvc[100]PDG particle code
        intNEipvert[300]PDG particle code X
->
        floatNEpvc[100][3]3-momentum (MeV/c)
        floatNEabspvert[300]Absolute momentum in the lab frame (MeV/c) X
->
floatNEpvc[100][3]3-momentum (MeV/c): consistent with fNeutrinoSim from ND::TTruthVerticesModule::TTruthVertex *vtxTrue
      //printf("test11 %d : %d %d %d \n", NeutrinoUtils::kPDGNumu,  rooVertex->StdHepPdg[0], rooVertex->NEipvc[0], rooVertex->NEipvert[0]);//test11 14 : 14 14 211 
      TVector3 tmpnu(rooVertex->NEpvc[0][0],rooVertex->NEpvc[0][1],rooVertex->NEpvc[0][2]);
      //printf("test12 %f %f\n", tmpnu.Mag(), rooVertex->NEabspvert[0]); //test12 1651.566620 -99999.898438 test2 1.651567; test12 935.936128 -99999.898438 test2 0.935936
      printf("test13 %f %f %f\n", tmpnu.Mag()/1e3, tmpnu.Phi()*TMath::RadToDeg(), tmpnu.Theta()*TMath::RadToDeg());//test13 2.474604 -113.921453 1.691795, test21 2.474604 -113.921453 1.691795; test13 0.546727 -120.282626 1.743435, test21 0.546727 -120.282626 1.743435
      */

    }
  }

  if(nfound>1){
    printf("AnaND280::GetNeutrinoParentPDG nfound>1!! %d\n", nfound); exit(1);
  }

  if(nfound==0){
    NeutrinoUtils::fNeutrinoParentPDG = NeutrinoUtils::kLost;
  }
  else{
    NeutrinoUtils::fNeutrinoParentPDG = parid;
  }
}

void AnaND280::SetMuonRec(const Int_t ib)
{
  Int_t ntot = 0;
  Int_t nfound = 0;

//  Int_t maxpTid = -999;
//  Double_t mpT = -1e10;

  Double_t mpp = -1e10;
  const ND::TGlobalReconModule::TGlobalPID * maxtrk = 0x0;

  // IN TRACK (PID) LOOP
  for(Int_t ii = 0; ii < fNPIDs; ii++){

    const ND::TGlobalReconModule::TGlobalPID *trTrack = (ND::TGlobalReconModule::TGlobalPID*)(*fPIDs)[ii];

    if(!BunchTime::IsPosTOK(trTrack->FrontPosition.T())){
      continue;
    }

    // #CUT#* BUNCH TIMING
    if(ib != BunchTime::GetBunchNo(trTrack->FrontPosition.T())){
      continue;
    }
    NeutrinoUtils::FillCount(fList, EventsPassedBunchTiming);
  
    // #CUT#* TPC INFO (if no TPC info for the track go to the next track (next j) in the entry
    if(trTrack->NTPCs == 0){
      continue;
    }
    NeutrinoUtils::FillCount(fList, EventsPassedTPCInfo);

    // #CUT#* VERTEX IN FGD1 FID VOL
    if(!FGDCut(&(trTrack->FrontPosition))){
      continue;
    } //if not in FGD1, reject
    NeutrinoUtils::FillCount(fList, EventsPassedVertexCut);
  
    // #CUT#* FORWARD-GOING TRACKS ONLY
    if(trTrack->isForward==0){
      continue;
    } //CUT forward tracks only
    NeutrinoUtils::FillCount(fList, EventsPassedForwardCut);

    Bool_t istpcok = kFALSE;
    Bool_t ismuon = kFALSE;
    IsTPCOK(trTrack, istpcok, ismuon);

    if(!istpcok){
      continue;
    }
    ntot++;
    NeutrinoUtils::FillCount(fList, EventsPassedTPC2Comp);

    if(ismuon){
      // Store momentum and position of highest momentum track
      if(trTrack->FrontMomentum > mpp){ 
        mpp = trTrack->FrontMomentum;
        maxtrk = trTrack;
      }
      nfound++;
      NeutrinoUtils::FillCount(fList, EventsPassedMuonCandCut);
    }
  }

  delete NeutrinoUtils::fMuonRec;
  NeutrinoUtils::fMuonRec = 0x0;

  ftTidPassMu = -999;
  ftRidPassMu = -999;

  //multi-muon allow                                                                                                                                                                                                    
  if(nfound>=1){
    TVector3 mup3(maxtrk->FrontDirection.X(), maxtrk->FrontDirection.Y(), maxtrk->FrontDirection.Z());
    mup3 *=  maxtrk->FrontMomentum * 1e-3 / mup3.Mag();

    const Double_t muE = TMath::Sqrt(NeutrinoUtils::MuonMass()*NeutrinoUtils::MuonMass()+mup3.Mag2());
    NeutrinoUtils::fMuonRec = new TLorentzVector(mup3, muE);

    ftTidPassMu = maxtrk->TrueParticle.Vertex.ID; //This is for the true vertex
    ftRidPassMu = maxtrk->TrueParticle.ID; //Trajectory ID for the track
  }

  NeutrinoUtils::SetMultiplicity(ntot);
  NeutrinoUtils::SetMuonNRec(nfound);
} 

Bool_t AnaND280::HasVetoTPCPOD(const Int_t ib)
{
  //Loop through the PIDs to see how many P0D, P0DECal and TPC1 objects there are in the bunch.
  //We want to veto events where any of these are non-zero.

  Int_t tpc1bb = 0, tpc2bb = 0, tpc3bb = 0, pdb = 0;

  // Loop over objects to describe a reconstructed PID
  for(Int_t ipid = 0; ipid < fNPIDs; ipid++){
    // Gets each PID object
    const ND::TGlobalReconModule::TGlobalPID *gbTrack = (ND::TGlobalReconModule::TGlobalPID*)(*fPIDs)[ipid];
    // Gets front position of track and checks if timing correct
    if(ib != BunchTime::GetBunchNo(gbTrack->FrontPosition.T())){
      continue;
    }
    // Checks number of POD objects for this PID object
    pdb = pdb+gbTrack->NP0Ds;
  }

  if(pdb){
    return kTRUE;
  }  
  NeutrinoUtils::FillCount(fList, EventsPassedPODVeto);


  for(Int_t itrk = 0; itrk < fNTracks; itrk++){
    //Loop over all tracks. Using trTrack to access Tracks[itrk] - each track object?
    const ND::TTrackerReconModule::TTrackerResult *trTrack = (ND::TTrackerReconModule::TTrackerResult*)(*fTracks)[itrk];
    
    if(!BunchTime::IsPosTOK(trTrack->Position.T())){
      continue;
    }
      
    // Feeding in track time to bunch calculator bunchtime.cxx
    // If bunch no. calculated from track position in time is not same as current bunch continue
    if(ib != BunchTime::GetBunchNo(trTrack->Position.T())){
      continue;
    }

    // Create TClones array TPC of TPC objects for each track?
    const TClonesArray * TPC = trTrack->TPC;
    // Loop over number of TPC objects per track
    for(Int_t itp = 0; itp < trTrack->NTPCs; itp++){
      // tpctrack accesses each TPC object and checks which detector it originated it
      ND::TTrackerReconModule::TTPCTrack *tpctrack = (ND::TTrackerReconModule::TTPCTrack*)(*TPC)[itp];
      if(tpctrack->Detector == 1){tpc1bb++;}
      if(tpctrack->Detector == 2){tpc2bb++;}
      if(tpctrack->Detector == 3){tpc3bb++;}
    }
  }
  
  if(tpc1bb){
    return kTRUE;
  }
  NeutrinoUtils::FillCount(fList, EventsPassedTPC1Veto);

  return kFALSE;
}

Int_t AnaND280::IniIO(const TString filelist)
{
  //Using bits of tracker and global recon - global generally better but tracker better for 1 bit...
  freconTracker     = NeutrinoUtils::InputFiles(filelist, "ReconDir/Tracker");
  freconGlobal      = NeutrinoUtils::InputFiles(filelist, "ReconDir/Global");           
  //fBasicHeader      = NeutrinoUtils::InputFiles(filelist, "HeaderDir/BasicHeader");
  //fBasicDataQuality = NeutrinoUtils::InputFiles(filelist, "HeaderDir/BasicDataQuality");
  //fBeamSummaryData  = NeutrinoUtils::InputFiles(filelist, "HeaderDir/BeamSummaryData");

  fmcTruthRooTrackerVtx = NeutrinoUtils::InputFiles(filelist, "TruthDir/NRooTrackerVtx");    
  fmcTruthVertices      = NeutrinoUtils::InputFiles(filelist, "TruthDir/Vertices");          
  fmcTruthTrajects      = NeutrinoUtils::InputFiles(filelist, "TruthDir/Trajectories");      

  //---
  //takes too much time! 
 // freconTracker->SetBranchStatus("*",0);
 // freconGlobal->SetBranchStatus("*",0);
 // fmcTruthRooTrackerVtx->SetBranchStatus("*",0);
 // fmcTruthVertices->SetBranchStatus("*",0);
 // fmcTruthTrajects->SetBranchStatus("*",0);
 //
 // freconGlobal->SetBranchStatus("NPIDs",1);
 // freconGlobal->SetBranchStatus("PIDs",1);
 // freconTracker->SetBranchStatus("NTracks",1);
 // freconTracker->SetBranchStatus("Tracks",1);
 // fmcTruthRooTrackerVtx->SetBranchStatus("NVtx",1);
 // fmcTruthRooTrackerVtx->SetBranchStatus("Vtx",1);
 // fmcTruthTrajects->SetBranchStatus("NTraj",1);
 // fmcTruthTrajects->SetBranchStatus("Trajectories",1);
 // fmcTruthVertices->SetBranchStatus("NVtx",1);
 // fmcTruthVertices->SetBranchStatus("Vertices",1);

  //---

  //fBasicHeader->SetBranchAddress("EventID", &fEventID);

  //NTracks matches the name in the oaAnalysis .root files and is the number of tracks in the tracker per MC entry
  //Int_t NVertices;
  fPIDs = new TClonesArray("ND::TGlobalReconModule::TGlobalPID", 200);
  //TClonesArray* GVertices = new TClonesArray("ND::TGlobalReconModule::TGlobalVertex", 200);
  //reconGlobal->SetBranchAddress("NVertices",&NVertices);
  //reconGlobal->SetBranchAddress("Vertices",&GVertices);
  freconGlobal->SetBranchAddress("NPIDs",&fNPIDs);
  freconGlobal->SetBranchAddress("PIDs",&fPIDs);
  
  fTracks = new TClonesArray("ND::TTrackerReconModule::TTrackerResult", 20);
  freconTracker->SetBranchAddress("NTracks", &fNTracks);
  freconTracker->SetBranchAddress("Tracks", &fTracks);
    
  fVtx = new TClonesArray("ND::NRooTrackerVtx", 100);  // GDB [GN]
  fmcTruthRooTrackerVtx->SetBranchAddress("NVtx",&fNVtx);   // GDB
  fmcTruthRooTrackerVtx->SetBranchAddress("Vtx",&fVtx);  // GDB
    
  fTrajectories = new TClonesArray("ND::TTruthTrajectoriesModule::TTruthTrajectory", 20);
  fmcTruthTrajects->SetBranchAddress("NTraj", &fNTraj);
  fmcTruthTrajects->SetBranchAddress("Trajectories", &fTrajectories);
    
  fVertices = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex", 20);
  fmcTruthVertices->SetBranchAddress("NVtx", &fNTVtx);
  fmcTruthVertices->SetBranchAddress("Vertices", &fVertices);

  //fBasicHeader->GetEntry(0);
  //BasicDataQuality->GetEntry(0);
  //BeamSummaryData->GetEntry(0);
  fmcTruthRooTrackerVtx->GetEntry(0);  // GDB
  fmcTruthVertices->GetEntry(0);  // GDB
  fmcTruthTrajects->GetEntry(0); //HMOK


  return freconTracker->GetEntriesFast();
}

void AnaND280::DoAnalysis(char *filelist, const Int_t ntarget, const TString tag)
{
  load();
  
  delete fList;
  fList =NeutrinoUtils::GetHistList();
  
  TTree *tout = NeutrinoUtils::GetTree();
  
  SetNeutrinoRec();
  
  std::cout<<"AnaND280::DoAnalysis ntarget "<<ntarget<<std::endl; 

  TStopwatch iotime;
  iotime.Start(kTRUE);

  const Int_t nentries = IniIO(filelist);  
  std::cout<<"filelist "<<filelist<<" "<<nentries<<std::endl;

  iotime.Stop();
  iotime.Print();
  
  Int_t nsel = 0;
  Int_t nfail = 0;
  Bool_t kstop = kFALSE;

  TStopwatch localtime;
  localtime.Start(kTRUE);

  for(Int_t jentry = 0; jentry < nentries; jentry++){                   //  *IN ENTRIES LOOP*
    const Int_t ientry = freconTracker->LoadTree(jentry);
    //Add event ticker
    if(ientry<0){
      printf("AnaND280::DoAnalysis ientry < 0 jentry %d\n", jentry);
      break;
    }

    if(jentry%10000==0){
      printf("%d/%d nsel %d nfail %d\n", jentry, nentries,nsel, nfail);
      localtime.Stop();
      localtime.Print();
      localtime.Start(kFALSE);
    }
        
    freconTracker->GetEntry(jentry);
    freconGlobal->GetEntry(jentry);
    //fBasicHeader->GetEntry(jentry);

    //fBasicDataQuality->GetEntry(jentry);
    //fBeamSummaryData->GetEntry(jentry);
    fmcTruthRooTrackerVtx->GetEntry(jentry);  // GDB
    fmcTruthVertices->GetEntry(jentry);  // GDB
    fmcTruthTrajects->GetEntry(jentry); //HMOK

    //--------------------------------------------------------------------------------------

    for(Int_t ib = 0; ib < 8; ib++){                       //  **IN BUNCH LOOP**
      NeutrinoUtils::FillCount(fList, EventsSeen);
      
      NeutrinoUtils::FillCount(fList, 0);

      //---------------------------------------------------------------------------------------------------------                                                                                                         
      //reconstruction                                                                                                                                                                                                    
      //---------------------------------------------------------------------------------------------------------                                                                                                         
      
      if(HasVetoTPCPOD(ib)){
        continue;
      }
      
      SetMuonRec(ib);
      
      NeutrinoUtils::FillCount(fList, 3);
      
      if(NeutrinoUtils::GetMuonNRec() == 0){
        nfail++;
        continue;
      }
      NeutrinoUtils::FillCount(fList, EventsPassedMuonCut);      
      NeutrinoUtils::FillCount(fList, 4);
      
      if(NeutrinoUtils::GetMuonNRec()>1){
        NeutrinoUtils::FillCount(fList, 5);
      }
      
      //std::cout << "Event: " << fEventID << std::endl;
      
      NeutrinoUtils::SetMuonPtRec();

      NeutrinoUtils::SetRecKinematics();

      //---------------------------------------------------------------------------------------------------------
      //simulation
      //as of r2685, the histogram output are all about simulation (key words: neutrinoE, Q2, sim)
      //if run on real data, no hist will be filled
      //tested by commenting out the following simulation part
      //---------------------------------------------------------------------------------------------------------

      SetNeutrinoParentPDG(ib);

      SetNeutrinoMuonSim();

      NeutrinoUtils::FillCount(fList, 1);

      NeutrinoUtils::FillCount(fList, 2);


      NeutrinoUtils::SetSimKinematics();

      tout->Fill();
      
      NeutrinoUtils::DoFillTH(fList);
      
      nsel++;
      
      if(ntarget>=0 && nsel>=ntarget){
        kstop = kTRUE;
        break;
      }
      
    }//end of bunch loop
    if(kstop){
      break;
    }
    
  }//end of entries loop

  localtime.Stop();
  localtime.Print();
  
  printf("nentries %d nsel %d nfail %d delta %d\n", nentries, nsel, nfail, nentries-nsel-nfail);
  
  TFile *fout=new TFile(Form("nd280AnaOutput_tag%s_%d_%d.root", tag.Data(), nsel, nfail),"recreate");
  tout->Write();
  fList->Write("lout",TObject::kSingleKey);
  fout->Save();
  fout->Close();
  delete fout;
  
}
